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Abstract 

After a brief comprehensive review of old and new results on the well 
known Fermi-Pasta-Ulam (FPU) conservative system of TV nonlinearly 
coupled oscillators, we present a compact linear mode representation of 
the Hamiltonian of the FPU system with quartic nonlinearity and peri- 
odic boundary conditions, with explicitly computed mode coupling coef- 
ficients. The core of the paper is the proof of the existence of one-mode 
and two-mode exact solutions, physically representing nonlinear standing 
and travelling waves of small wavelength whose explicit lattice represen- 
tations are obtained, and which are valid also as N ^ oo. Moreover, and 
more generally, we show the presence of multi-mode invariant subman- 
ifolds. Destabilization of these solutions by a parametric perturbation 
mechanism leads to the establishment of chaotic in time mode interac- 
tion channels, corresponding to the formation in phase space of bounded 
stochastic layers on submanifolds. The full mode-space stability problem 
of the N/2 zone-boundary mode is solved, showing that this mode be- 
comes unstable through a mechanism of the modulational Benjamin-Feir 
type. In the thermodynamic limit the mode is always unstable but with 
instability growth rate linearly vanishing with energy density. The phys- 
ical significance of these solutions and of their stability properties, with 
respect to the previously much more studied equipartition problem for 
long wavelength initial excitations, is briefly discussed. 
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Keywords: Anharmonic lattices; Energy equipartition; Periodic solutions 
Hamiltonian bifurcations; Stochastic layers. 
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1 Introduction 

The numerical experiment by Fermi, Pasta, Ulam (FPU) and Tsingou|^ 
perhaps the first of aU times, was the first attempt to check the predictions of 
classical statistical mechanics concerning the dynamics of a Hamiltonian system 
of coupled oscillators with a large number N of degrees of freedom. The result 
of this experiment was a big surprise: the expected relaxation to equipartition 
of energy among the linear normal modes was not revealed, during the time 
of observation and with low energy initial excitations. This also implies that 
ergodicity is not an obvious consequence of the non-existence of analytic first 
integrals of the motion besides the total energy (and possibly total momentum). 
After an initial growth of the energy in the neighbouring modes (as expected) 
they observed that the energy sharing was restricted only to the first modes, 
with a quite regular dynamics, rather than detecting a gradual and continuous 
energy flow from the first excited mode to the higher ones and a stochastic 
dynamics. Even more surprisingly, at later times, almost all the energy was 
flowing back into the initially excited mode, so that the system seemed to possess 
quasiperiodicity properties. Chirikov et al. |^ showed later, also with numerical 
experiments that, at sufhciently high energies, the FPU model did relax 
to the equipartition state, on times which become smaller and smaller as the 
energy is increased. It then became clear that such a system was displaying 
qualitatively different behaviours as the energy, fixed by the initial condition, 
was varied. Chirikov and Izrailev Q gave also the first linear normal modes 
representation of the system and applied the nonlinear resonances theory to 
predict the existence of the transition between the two behaviours. They did also 
recognize the resemblance between these results and the Kolmogorov-Arnol'd- 
Moser theorem, first formulated, although not completely demonstrated, by 
Kolmogorov just in the same year of the FPU numerical experiment. These 
results successively stimulated a lot of numerical studies aiming at determining 
the dependence of the different observed behaviours of the FPU system on the 
number N of degrees of freedom (see Ref. ^ and references cited therein, and 
also 10 I, I). 

It is now well established that the transition between a quasi-integrable 
behaviour to a mixing one, in the FPU system as well as in similar models 
(Lennard- Jones in 1-D and 2-D, (p^), is controlled by the energy. At small ener- 
gies the motion is weakly chaotic with positive but small Lyapunov exponents, 
revealing the presence of thin stochastic layers in the phase space, which is 
mostly filled with KAM tori. At higher energies the maximum Lyapunov expo- 
nent and the Kolmogorov- Sinai entropy rise considerably, revealing the growth 
of stochastic regions. In this energy range the temporal signals of normal mode 
energies also become chaotic. Well above the transition region all the signatures 
of large-scale chaos are present, with TV — 1 {N — 2 for systems with periodic 

^M. Tsingou contributed to the numerical work and then did not participate in the writing 
of the report. 
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boundary conditions) positive Lyapunov exponents, fast diffusion of the orbit 
over the constant energy surface, rapidly decaying spatio-temporal correlation 
functions (also between the modal energies). It has been also shown that, in 
this region, a density function of Lyapunov exponents can be introduced in the 
thermodynamic limit, which can also be analytically estimated with a random 
matrix approximation |Tl|] . 

On the other hand, the quasiperiodic behaviour at low energy was interpreted 
in terms of a continuum limit of the model. It was in fact observed by Zabusky 
and Kruskal [^2] that, increasing the number N of oscillators while holding 
constant the length of the chain, the FPU model with cubic nonlinearity can be 
reduced to the Korteweg-deVries (KdV) equation in the limit of low amplitude 
and long wavelength excitations, after selecting one of the two possible directions 
of propagation for travelling waves. For the KdV equation the same authors 
found numerically the existence of localized solutions which preserve their shape 
even after mutual collisions. To emphasize the extremely stable and localized 
character of such solutions Zabusky and Kruskal called them solitons. It was 
shown later ||l^, |l^ that the KdV equation is a completely integrable infinite- 
dimensional Hamiltonian system, whose initial value problem can be solved by 
the spectral tranform method [|6[ |l^ . The quasiperiodic behaviour of the FPU 
system at low energies is consequently believed to be the result of the presence on 
the lattice of several solitons emerging from the long wavelength initial condition 
and travelling with different speeds so that, at particular times, they reproduce 
almost exactly the initial condition. Quasiperiodic behaviour is destroyed at 
larger energies, above the threshold of transition to equipartition A deeper 
understanding of this transition phenomenon requires a systematic study of its 
dependence on the number of degrees of freedom, on the initial condition and 
on the time of observation. Such a study has not been fully completed yet. 

1.1 The FPU model 

The FPU system is an approximate model for the behaviour of a classical solid 
at low temperatures. The reduction of complexity in comparison to the real 
physical situation is considerable. First of all only one spatial dimension is 
considered and then the interaction (typically of the Lennard- Jones kind) is 
expanded in the small displacements around the equilibrium positions of the 
atoms or molecules, i.e. the weakly anharmonic case is considered (in prac- 
tice the case of rather low temperatures). As a consequence, several problems 
must be discussed concerning the physical relevance of the FPU model. We 
must first ask if the observed phenomena are peculiar of the model or can be 
instead revealed in a wider class of more physical models. This issue is now 
clarified and we can state with no doubt that different models in one and two 
dimensions share the same phenomenology (see Ref. |^ for a review). In par- 
ticular, practically all the models which, as the FPU one, are constituted by 
weakly coupled anharmonic oscillators, show the prevalence of ordered motion 
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at sufficiently low energies, while at higher energies their phase space becomes 
chaotic. An important point is that the transition between the "ordered" re- 
gion and the "chaotic" one happens in an energy range relevant for physics. 
For example, in a model of oscillators in two dimensions with Lennard-Jones 
interactions , the transition, with the physical parameters of Argon, is found 
to happen around a temperature of 5 °K. Some recent molecular dynamics cal- 
culations |l9|| on a model for a Xenon crystal with diluted impurities made of 
Iodine molecules, show that the time required for the equipartition of energy 
between the initially excited I2 molecule and the Xe matrix abruptly increases 
at approximately 30 °K; such a physical situation could be investigated in a 
real experimental sample. 

Having established the physical relevance of the phenomena studied in the 
FPU model, we proceed to its definition. The FPU model is a one-dimensional 
chain of oscillators with unit mass, with weakly nonlinear nearest-neighbour 
interaction (the lattice spacing is also taken of unitary length). Calling and 
Pi the coordinates and, respectively, the momenta of the oscillators, the model 
is defined by the following Hamiltonian: 



N ■) N 



2 r 



(1) 



where qi = qn+i and r — 3, for the so-called FPU-a model, while r = 4 for 
the FPU-/? model. We consider in the present work only the FPU-/3 model, 
with periodic boundary conditions to allow travelling wave solutions (whereas 
in the original FPU paper fixed ends boundary conditions were chosen) . At fixed 
energy, the coupling constant (3 determines the amount of nonlinearity in the 
model. Conversely, as it is more natural, for a fixed f3 the increasing departure 
from the harmonic behaviour is controlled increasing the energy. It can be easily 
shown that the dynamics depend only on the parameter f3E'^^/'^~^\ where E is 
the total energy, fixed by the initial condition. 

Hamiltonian (|l|), written in linear normal coordinates {Qk,Pk) (phonons) 
becomes 

H=\Y.(.Pk+^lQl)+f3ViQ) . (2) 

k 

with frequencies tUk = 2sin(7rfe/A^) in the case of periodic boundary conditions 
(the harmonic frequency spectrum is obvioulsy different from the case of fixed 
boundary conditions studied by FPU) and we have LUk = uJN-kj so that there 
are only N/2 different frequencies (if N is even, for simplicity). The harmonic 
energy of mode k is defined by Ek ~ (-P| + ll}IQI)/2. 

The FPU experiment aimed at showing the progressive decorrelation of the 
system during its temporal evolution, eventually leading to ergodic behaviour. 
To this end they chose an initial condition far from equilibrium, giving all energy 
to the lowest (k — 1) normal mode only, and then calculating the instantaneous 
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energies Ek{t) of all modes. They expected to see a progressively uniform redis- 
tribution of energy between all modes, caused by the small anharmonic coupling 
between them. On the contrary they observed the well known FPU recurrent 
behaviour: energy was flowing back regularly to mode k = 1 after an initial 
share. Return to the initial condition is not exact, but the possibility that re- 
laxation is present on longer times was ruled out by the numerical experiment 
of Tuck and Menzel |Q , who first observed the "superperiod" phenomenon. 

The behaviour of average energies Ek{T) = Ek{t)dt versus T, shows 

that the systems relaxes to an asympotic state, but the latter is different from 
that in which equipartition of energy holds. On the contrary, at higher en- 
ergy the equipartition state is reached in a relatively short time. A transi- 
tion is present from an "ordered" state, in which the system seems not to be 
relaxing towards equipartition, showing recurrent behaviour in time, to a differ- 
ent one where, on the contrary, equipartition is quickly reached, which we call a 
"chaotic" state (though these terms are improper because also in the "ordered" 
state there are chaotic solution, but they fill a part of small measure in the 
phase space). 

Various indicators have been employed to characterize this transition, either 
based on the distribution of energy among the normal modes , or on the local 
divergence rate of nearby trajectories in phase space (Lyapunov exponents) [Q. 
In particular, the first clear numerical evidence of the transition to equiparti- 
tion (or "chaotic") state was obtained in Ref. Q by the use of an appropriate 
Shannon entropy S defined in the space of modes ("spectral entropy"). At a 
given value of the control parameter 

where E is the total energy, the "spectral entropy" was shown to increase, 
reaching then its maximal value. 

Despite these results the problem is still far from beeing understood, es- 
pecially concerning the dependence of the phenomenology on the number of 
degrees of freedom and on the time scales of observation. 



1.2 Recent results 

We will focus on the results which are valid at large N (thermodynamic limit). 
The first result we want to mention concerns the Lyapunov spectrum well above 
the region of transition to equipartition (e ^ 1). There are strong numerical 
evidences supporting the possibility to define a function which repre- 
sents the density of Lyapunov exponents in the thermodynamic limit, in the case 
when the only vanishing exponents are those related to space and time trans- 
lational symmetries (momentum and energy conservation, respectively). After 
labelling the Lyapunov exponents in the decreasing order Ai > A2 > . . . > A2JV, 
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the density function (j){x) is given by 



0(x) = lim X^n{N) . (4) 

The existence of this hmit was conjectured by RueUe for the Navier-Stokes 
equation, and verified by ManneviUe for the Kuramoto-Sivashinsky equa- 
tion. The existence of such a function makes it possible to define (using the Pesin 
relation) the Kolmogorov-Sinai entropy as an intensive quantity in the thermo- 
dynamic limit, thus proving the persistence of a state of large-scale chaos also 
when the number of degrees of freedom goes to infinity. So, when the intensive 
control parameter e is big enough, the trajectory of the system visits in a short 
time all accessible regions of phase space and, as a consequence, equipartition 
of energy among the normal modes is quickly reached. The evidence for the 
existence of the Lyapunov exponents density function shows, in an indirect way, 
that equipartition of energy in the FPU model occurs also in the thermody- 
namic limit, when the energy density is high enough. Unfortunately there are 
no analytical results proving the existence of the function 0(a;), apart from those 
obtained by Eckmann and Wayne using a random matrix approximation [|ll|. 
A recent discussion of this problem was also given by Sinai , who has indeed 
obtained a rigorous proof of the existence of an analogous quantity, obtained 
interchaning the order of the limits TV — > oo and t —> oo (the latter appearing 
in the definition of the Lyapunov exponents). Though this is not, a priori, the 
same thing as (j), it can be argued to be more relevant for statistical me- 
chanics; moreover it has a closer correspondence to the numerical procedures 
used to determine (j). 

Pettini and Landolfi Q found a rapid decrease of the maximum Lyapunov 
exponent Ai as the control parameter e is decreased below a value e ~ 0.1; such 
a behaviour does not change significantly when N is increased. Even if this 
does not necessarily prove the onset of large-scale chaos (because only a few 
directions, among many, might be responsible for the increase of the divergence 
rate), the persistence of the result when N is varied reveals the presence of an 
important phenomenon occurring in phase space. This result is consistent with 
the results obtained using the "spectral entropy" , which also sharply increases 
around e « 0.1 (although a slower diffusion of the entropy with time is also 
present for smaller values of e) |2^, . 

More recently it has been suggested that the largest Lyapunov exponent 
can be estimated from the temporal evolution of the vector ^ representing the 
deviation between geodesies on the Riemannian manifold constructed with the 
Eisenhart metric over the enlarged configuration space-time. In particular, it 
was shown that the norm ■(/> = ||^|| obeys the equation 

4, = -Kr^P (5) 

where Kr — V'^V{q)/N is the positive Ricci curvature of the manifold and 
V the potential of the FPU model. Supposing that the quantity Kr and the 
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magnitude of its fluctuations can be calculated, in the thermodynamic limit, 
using the Gibbs equilibrium measure and solving for the instability growth rate 
of Eq. (|^), Casetti, Livi and Pcttini [ p8[ obtain values of the maximum Lyapunov 
exponent compatible with those reported in The exponent would be zero 
in the absence of temporal fluctuations of Kn. It is an intriguing fact that the 
time-periodic dependence of Kji on stable periodic (elliptic) orbits must differ 
enough from that on unstable periodic (hyperbolic) orbits, even if they are close 
in phase space, to make the Lyapunov exponent vanish. In fact, the positivity 
of Lyapunov exponent is caused, in this case, by the fact of being out of the 
spectrum of the periodic operator in (||). A more complicated case is the one of 
"chaotic" functions Kr. It is treated in Ref. and has a close resemblance 
to the Anderson localization problem in one dimension. 

Concerning perturbation theory results, Galgani, Giorgilli and collaborators 
in a recent paper |9| summarize the findings obtained using Nekhoroshev esti- 
mates The latter permit to evaluate upper boimds for the time variation 
of the unperturbed actions on times that, though being finite, increase expo- 
nentially as the perturbation parameter is decreased (in the FPU case the /3 
parameter of the nonlinear term). It is possible, using this approach, to find 
results valid for initial conditions on open sets in the phase space, as opposed to 
methods based on the Kolmogorov-Arnold-Moser theorem (on the other hand 
the latter has the advantage to give statements valid for all times) . The stability 
time r of the single unperturbed actions (or action "freezing" time) is found to 
scale as 

r = r.exp(^)^ (6) 

where, in general, both t, , /3* and d depend on N. The most important de- 
pendence on N is that of d: the best estimated so far obtained for FPU gives 
d ~ N~^, a result confirmed also by numerical simulations (so the estimate 
seems to be optimal). This result suggests that in the thermodynamic limit 
the freezing times might become short, or even vanishing, and the region of 
violation of energy equipartition could disappear. We must, however, remem- 
ber that such estimates are valid in an energy region shrinking to zero as N is 
increased, and moreover they are lower bounds, so they could be irrelevant. A 
positive result was however obtained on a modified FPU model with alternating 
light and heavy masses (like in a diatomic solid). In this case the harmonic 
frequency spectrum splits in two separated components: an "acoustic branch" 
with low frequencies and an "optical branch" at high frequencies; the latter 
is almost completely degenerate as the ratio between the masses is increased 
and N ^ oo. Then it is possible to give Nekhoroshev estimates for subsets of 
unperturbed actions as the total harmonic energy of the optical modes, which 
turns out to be frozen on a time given by the law (^, but now with d — 1, 
r* ~ and /3* ~ N^^. Such analytical estimates suggest again the vanish- 
ing of the freezing time in the thermodynamic limit, but in this case numerical 
simulations show that the estimates are far from being optimal: seems to be 
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independent of N, while r* has a weak dependence on N, possibly like (In A^)"^. 
The authors of Ref. then argue in support of the relevance of Nekhoroshev's 
like estimates in the thermodynamic limit. Unfortunately, for the FPU model 
it is not easy to identify a natural decomposition in subsets of unperturbed ac- 
tions. Perhaps this might be possible in the light of the present work, where we 
show the existence of subsets of normal modes where energy remains trapped for 
suitable initial conditions, so defining reduced (integrable or not) Hamiltonians 
on submanifolds of the constant energy hypersurface. It would be interesting 
to attempt Nekhoroshev-like estimates when the energy is restricted to these 
submanifolds. 

More recently De Luca, Lichtenberg and Lieberman studied the ap- 
proach to equipartition in the FPU model exploiting normal form theory to find 
an effective Hamiltonian describing the interaction among a reduced number of 
modes (four in their treatment). The main result of their theory is that above 
a critical energy Ec the system reaches a near-equipartition state, in a time 
proportional to iV^; below this critical energy the time needed increases even 
faster with N (perhaps exponentially). This holds when the initial excitation 
is given to a subset of low modes whose number does not increase with N. If 
instead the excited modes are a subset scaling proportionally to N {k (x N), 
the typical time scale to quasi-equipartition increases like N. These predic- 
tions are also supported by numerical simulations. The idea to construct the 
reduced Hamiltonian is to make a large N expansion of the dispertion relation 
LOk = 2sin[7rfc/2(iV + 1)] ~ [Trk/{N + 1)] + 0{{k/N)^) (we are treating now 
the system with fixed ends) and to identify the four-wave resonance relations 
(fci-|-A;2-l-fc3-l-A;4 = 0) producing in the resonant normal form those angles which 
are slowly (adiabatically) varying; these latter are found to be dg — O1 + 63 — 26*2 
and 9sp = 62 + 64 — 26*3. The parameter controlling the deformation of the 
actions (monotonic in the energy) is seen to be proportional to NE, while the 
angles and dgp are slowly evolving with the frequency fl ^ kE/N'^; the lat- 
ter determines the characteristic evolution time for the actions r ~ N"^ /E for 
k ~ const., while r ~ N/E ii k (x N. The energy transfer to higher modes is 
present, but takes place on much longer times. Actually, it is also known that 
the energy fraction transferred to the highest modes is exponentially small in 
mode number 0. De Luca et al.'s idea is to use resonant normal forms 
and expand the dispersion relation for large N to obtain an effective Hamilto- 
nian. This latter expansion was also suggested by Shepelyansky however 
his purpose was to estimate the largest Lyapunov exponent which he was able 
to show, by his method, to be positive even at very small energies. So, accord- 
ing to Shepelyansky, the critical energy for the onset of chaos goes to zero as 
N is increased (in agreement with what is obtained in |Q). This last result is 
consistent with the numerical result by Pettini and Landolfi but the scaling 
law suggested by Shepelyansky is not consistent with numerical results. 

It is still difficult to arrange all these results in a single coherent framework, 
and especially it is not clear if it can be really justified to neglect the interaction 
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with high modes. In fact, we show in this paper that the structure of the 
interaction among the modes is rather complex. 

Kantz, Livi and Ruffo ||2^, indipendently from De Luca et al., revealed 
with numerical experiments that the typical time scale for equipartition in- 
creases with TV. More precisely, considering the functional dependence of the 
"spectral entropy" S on energy E, time t and on the number N of oscilla- 
tors, these authors found that S depends only on the two variables E and t/N: 
S{E,t,N) — f{E,t/N), for initial conditions with k oc N. Moreover, for ini- 
tial conditions with k ~ const, they found S{E,t,N) = g{E/N,t/N). This 
last result is in apparent contradiction with De Luca et al., who rather seem 
to suggest S{E,t,N) ~ g{E,t/N'^): but if both results have to be valid, it is 
clear that S must be a function of the single parameter Et/N"^ . This prediction 
has been recently numerically confirmed ||2^ . Moreover, a completely indepen- 
dent analytical method based on the derivation of the breakdown (shock) time 
for the non-dispersive limit of the mKdV equation (to which the FPU-/3 model 
reduces in the continuum limit) leads to the same prediction [ |33t . However, a 
correction to this scaling theory due to a phase-space filling diffusive process 
is also numerically revealed pa, which finally gives a dependence of S on the 
single parameter j3Ekt/ (N'^wN). It should be observed that even in the ther- 
modynamic limit at fixed E/N and with k oc N, will the time scale of the 
system then diverge as y/N. This is very intriguing because it suggests (also in 
the spirit of Ref. |p3[) that far from equipartition states might persist in time 
if we perform the N oo limit before the t ~f oo limit (a quite uncommon 
procedure for equilibrium statistical mechanics, but perhaps a meaningful one 
in non-equilibrium). 

The purpose of this paper is twofold. First, we want to show several exact 
periodic and quasiperiodic solutions for the FPU-/3 system with periodic bound- 
ary conditions, whose majority appears to be new. Those which have already 
appeared in the literature are rederived by a unified method, which consists es- 
sentially in "extending" some harmonic modes to the anharmonic case. Second, 
having at our disposal explicitly known solutions, we take a first step in the 
study of their stability and its connection with the onset of chaos in the system. 
To achieve these results we rely on the structure of modal interactions in the 
FPU-/3 model, which is derived in Section In Section || we discuss the specific 
fully integrable case of the three particle FPU-/3 model. In Section |4| we derive 
explicit one-mode and two-mode solutions. In Section ^ we identify two-mode 
and three-mode invariant submanifolds where the motion displays transitions to 
chaos after successive bifurcations. Section [s] is moreover devoted to a study of 
the restricted stability in a two-mode submanifold and to the full phase-space 
stability analysis of the k = N/2 zone-boundary solution. 
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2 Linear normal coordinates 
for the FPU-/5 system 

Hamiltonian (|l|) of the FPU-/3 system with N particles and periodic boundary 
conditions can be splitted in two contributions: 

H = Ho + Hi (7) 

with 

1 ^ 1 ^ 

^0 = 9 II Pn + 9 I](9n+1 - Qnf , (8) 



2^"" 2 

71—1 n— 1 



^1 = ^5](g„+l-g„)^ (9) 

where ^at+i = qi and /3 > 0. 

While the interaction in physical space has the simple structure of nearest- 
neighbour coupling, it is clear that it becomes highly tangled when it is rep- 
resented in modal (i.e. Fourier) space. Explicit expressions for the equations 
of motion in normal coordinates appeared sporadically in the literature for the 
FPU-/3 system in the fixed boundary case, and they were used for varied pur- 
poses. In Ref. Ql the modal equations were treated with the Krylov-Bogoliubov- 
Mitropolsky averaging technique, and then subjected to the criterion of overlap- 
ping resonances |Q to obtain estimates for the threshold ofstochasticity. Re- 
cently, in ^ the modal equations are the starting point to develop approxi- 
mations, based on normal form theory, to study the time scales to cquipartition 
in the case when the energy is initially in a single or small group of low-frequency 
modes. Estimates for the FPU recurrence period were obtained by perturbative 
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treatment of the modal equations 1 35 . Recently, ShoU and Henry applied 
a shifted frequency perturbation scheme to the modal equations to perform 
calculations of superperiod recurrence times, at low anharmonicity, in FPU-a 
and FPU-/3 chains with fixed ends. In Ref. the case of the FPU-/? chain 
with fixed ends and 15 moving particles excited from rest in the 11th mode 
was considered (this is the most studied case were the so-called induction period 
phenomenon |39j shows up); a detailed description of the early phases of 
the dynamics in the chain, resulting from numerical integration, was given and 
interpreted by means of the structure of modal couplings. 

In our treatment we need to use some selection rules arising from the struc- 
ture of the equations of motion in Fourier space. The first description of selection 
rules for the FPU-/? chain with fixed ends was given by Bivins, Metropolis and 
Pasta |39[ and later in more detail and for more general interaction potentials 
by ShoU [|0[ |l|. The case of periodic boundary conditions is slightly more 
involved, due to the degeneracy of the harmonic frequency spectrum, and to 
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our knowledge it has been given only a very brief treatment, with some missing 
points, in unpublished work by Sholl It is not however our purpose to 

give here a complete description of the selection rules arising in the periodic 



boundary case, which could be given along the lines of Refs. |^0|, 41 , and we 
will introduce only what is strictly necessary to derive our results. 



2.1 Normal coordinates for the periodic chain 

Although the transformation to normal coordinates of the harmonic Hamilto- 
nian Hq is quite an elementary topic, it is useful to consider it in some detail to 
fix the notations and make some remarks. 

The "unperturbed" Hamiltonian Hq can be written, introducing the column 
vectors q and p of coordinates and conjugate momenta, as 

Ho = ^'pp+^'qAq (10) 

where *• denotes transposition and A is the symmetric N x N matrix with 
elements 

i = 1,. ..,N 

Aij = —6i-ij + 2Sij — Si+ij — Si^i6j,N — Si^NSj^i , j — 1 ' N ^^^^ 

Let S be a iV X iV real orthogonal (*SS = I) matrix which makes A diagonal 
through the similarity transformation A i-^ *SAS. Then we can perform a 
canonical transformation from the old set of coordinates 

( *g; *p) = (91, ■ • ■,qN;pi, ■ ■ ■ ,pn) 

to the new set of coordinates 

( *Q; *P) = iQo,...,QN-i;Po,---,PN-i) , 
by the generating function 

N-l N 

F(9,P)= ^^5„fcPfcg„ . (12) 

k=0 n=l 

We have shifted the subscript of the new coordinates for reasons of convenience, 
i.e. to label the center of mass motion by the zero index. Through this trans- 
formation the harmonic part Hq takes the form of the Hamiltonian of a system 
composed of A^ — 1 uncoupled harmonic oscillators (the harmonic normal modes) 
plus a free particle (representing center of mass motion) . 

In view of the application of this transformation to the full anharmonic sys- 
tem it is important to remark that in the present case of periodic boundary 
conditions, owing to the degeneracy of the spectrum of A, the matrix S can be 
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chosen in infinite different ways, leading to different sets of normal coordinates 
which, at variance with the harmonic case, are not equivalent for the full anhar- 
monic Hamiltonian. Therefore the set of normal coordinates must be carefully 
specified. 

Defining, for integer k, 

a;fc = 2sin(^^^ , (13) 

the eigenvalues of A are /ife = ojf, for k E {0, 1, ... , [A^/2]|0, and they have 
multiplicity two, apart from /iq (=0) and (for even N only|^ ) which are 

non-degenerate. 

The uniquely defined (apart from a sign) normalized eigenvectors u*-*^) and 
u(^/2) corresponding respectively to /iq and have components 

ul^^ = n=l,...,iV (14) 

uW^) = i-^, n = l,...,A^ (15) 



Because of the degeneracy of the spectrum, there is freedom in the choice of 
an orthonormal basis in each two-dimensional eigenspace. Defining for k = 
1, . . . , [{N — l)/2] the vectors u^'^^(7), containing the arbitrary real parameter 
7, as 

^'nHl) = \flsu.(^-^+^ , n = l,...,N (16) 



it is easy to verify that 

Auf'^) (7) = ^IkU^''^ (7) , V7 e i2 (17) 

and that the set of N vectors {'w^''\'^)}ke{o,...,N-i} given by 

{u(0) for fc = 

uW(7) forfc = l,...,[i^] 

uW2) for k = N/2 ^^^> 

u(^-fe)(7 + f) for A:= [f] + l,...,iV-l 

is an orthonormal basis composed by eigenvectors of A (for each arbitrarily fixed 
7). The eigenvectors corresponding to the eigenvalue fit (k = I, ■ ■ ■ , [(A^— 1)/2]) 

^We denote by [x] the integer part of x, so that [A'^/2] = N/2 if A'' is even while [N/2] = 
(N - l)/2 if N is odd. 

^To keep generality while avoiding frequent distinctions, from now on we keep the conven- 
tion that every proposition where the value k = N/2 appears is intended to be valid only for 
the even TV case, while it can be simply omitted for odd A'^. This is true, in particular for (hS) 
and @. 
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are w(''')(7) = u('^^(7) and w^^ ^^^7) = u(''')(7 + §)• Once the basis, i.e. 7, is 
chosen, the matrix S^"*"^ which diagonahzes A is given by 



Using this matrix in the generating function of the form (|12|) , the set of normal 
coordinates and momenta {Q^'^\ -P*'''''), associated to the chosen basis, is defined 
by the canonical linear orthogonal transformation 

For each fixed k, the set of components {wl^^ {7)}ne{i,....N} 1 viewed as a fmiction 
of the lattice position n, gives the spatial pattern associated to the presence of 
the normal coordinate Q^^\ since = J2k=QQk'\^)^"'^i"f)- Choosing a 

basis amounts to a choice of the spatial phase factor 7 for the stationary basis 
waves used to Fourier analyze the spatial pattern of the chain, given by the 
set {qn}- Different basis give different Fourier components, i.e. different sets 
of normal coordinates. The coordinate Qq is proportional to the displacement 
of the system's center of mass, while for k g {1, . . . , [^^^j^]} both coordinates 

(Qi7\QiV-fc) correspond, in the Fourier analysis of the spatial pattern, to the 
presence of a component of wavelength N/ k (the lattice spacing being the unit 
of length) but they refer respectively to a sine and cosine wave (with respect 
to the chosen spatial phase). Finally the Qn/2 coordinate corresponds to the 
pattern of wavelength two, where adjacent particles have equal but opposite 
displacements. 

Every choice of the basis is equivalent as far as the linear system is concerned 
since Hamiltonian (|^) takes the same form in every set of normal coordinates 
(Q(7)^p(7)), This is no longer true for the full nonlinear Hamiltonian H, when 
the same canonical transformation is applied. 



2.2 Transformation of the full Hamiltonian 

We now choose once for all the basis with 7 — 7t/4. This is done for two reasons. 
Firstly, this choice simplifies the calculations since all the elements of the matrix 
S = S^^/"*-* are given by the single expression 

n = l,...,N 
k = 0,...,N -I 

(21) 

Moreover, we shall see that some of the normal modes of this basis can be 
"extended" to the anharmonic case. From now on we denote simply (Q, P) = 



1 




f 2'Kkn\ 


( 2iikn\ 




sin 


K N J 
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The quartic term Hi in is transformed into 

N /N-1 \4 
n=l \fc=0 / 



with SN+i,k = 5'i,fc- 
Since 



_ 2 fT:k{2n + l) ^ tt 

•Jn+l^fe ~ ^nk — \ Tr^k COS 



AT V 4 

it is useful to introduce the quantities 



f/4 ^ cos ( ^-j^ >- + 6 



which obey the following rules 



2 



The nonlinear term (^2|) can be rewritten in terms of these quantities 

N~l N 
i,j,k,l — l n—1 

Using the properties (^,(|2^) we get 

N 



n=l 

g / , V ^ n,i+j+k+l ^ ^ n,i+j-k-l ^ ^ n,i-j+k-l ^ 
n—1 

\ n,i+j+k-l ^ '-^n,i+j~k+l ^ ^ 1 



/2 -Ifl"^ 
-\-k — l ^ ^ ~k-\-l ^ ^ n,i—j+k+l n,i—j~k- 



n=l 

Using the identity, valid for integer r 

N 



n=l 

we obtain 



/. 7rr(27i + 1) \ _ / (-l)™iV if r = miV with m e Z 
exp yi ^ / 1 otherwise , 



N 
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and 



being 



N 



Y^Kr^NAr ,VreZ, (31) 



n=l 



^ _ J (—1)™ for r = mN with m E Z . , 

^ I otherwise . ^ ' 



From (M), im and (Bl^ we get 



with 



Hi{Q) — ^ LOitOjUJkUJiCijkiQiQjQkQi , (33) 

i^j,k,l—l 



Cijki — — Ai+j+fc+i + Ai+j-k-i + Ai-j+k-i + Ai^j^k+i (34) 



These integer coefficients are invariant under any permutation of the indices. 

The full Hamiltonian of the FPU /3-model for N oscillators with periodic 
boundary conditions in the new coordinates is 



1 1 

H{Q, P) - -P2 + _J2iP'+ ^iQi) + iii{Q) ■ (35) 



Eliminating center of mass motion, we get easily also the equations of motion 
for the remaining iV — 1 degrees of freedom, which, in second order form, read 

Qr = FriQl,..., Qn-i) , r^l,...,N-l 

FriQl, ... ,Qn-i) = ^^rQr — ^^J^ ^ UjUJkLOlCrjklQjQkQl , ^ ^ 

j,k,l=l 

where Fr is the generalized force in normal coordinate space. 



3 Intermezzo: the integrable = 3 case 

It is well known that the FPU /3-model with = 3 and periodic boundary 
conditions is integrable . This is due to the presence of a third independent 
integral of motion besides energy and total momentum. In ref. it is in fact 
shown that the quantity 

Miq,p) =pi{q2 - qs) +P2(93 - Qi) + Psiqi - 92) (37) 

is a constant of the motion. Our canonical form of the FPU Hamiltonian ( |33| , ^5| ) 
allows a simple interpretation of this new integral and clarifies why it is typical 
of TV = 3. 
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In the case A'^ = 3 there are only two non-zero harmonic frequencies, the 
coincident pair ujx = u}2 = \/3- Invariance under permutation of indices of 
the CijM coefficients in formula (|3^) allows to restrict the search for non zero 
coefficients among those with i > j > k > I. These are Ci^i^i^ = C2.2,2.2 = 3 
and C2,2.i,i = 1- The resulting Hamiltonian is 

H{Q,P)='-P^ + l{P? + P^) + l{Ql + Ql) + l(3{Ql + Qiy . (38) 

Obviously the center of mass motion is free, and Pq is conserved. More im- 
portant is the invariance of the potential term in ( |3^ ) under rotations in the 
(Qi, Q2) plane, which is transparent in these new coordinates. This means that 
the pseudo-angular momentum L — {Q1P2 — Q2P1) is conserved. Going back 
to old coordinates and momenta by the inverse of the transformation (pO|), we 
get easily the new integral in (|37|): 

M = V3L . (39) 

Thus, the presence of this new integral is not a surprise since it arises from 
a simple symmetry which was not apparent in the original coordinates; this 
symmetry is no more present for > 3, since the potential does not have in 
general rotational invariance in the planes (Qfc, Qn-r), but it holds true in the 



invariant subspace {N/3,2N/3} (see Subsection 4.2). 

We finally observe that one can give a different interpretation of the sys- 
tem considered in this section. Let us think of the coordinates (<Zi,(72,93) as 
cartesian coordinates of a single particle governed by the Hamiltonian (^) with 

= 3, thus identifying the configuration space of the three one-dimensional 
oscillators system with that of a single particle in three-dimensional space. The 
transformation to normal coordinates, being induced by an orthogonal matrix. 



is a rotation of the reference frame in this context, and the form (38) taken 
by the hamiltonian in these new coordinates reveals the cylindrical symme- 
try of the potential around the w*^"^ axis and its independence from the Qo 
coordinate associated to the latter. From this immediately follows the conser- 
vation law for L, the angular momentum component along the w^^^ axis, and 
for Pq = {pi +P2 +P3)/VS which is now interpreted as the particle's momentum 
component p ■ w*^"^ ; in fact there is no force along that direction. 

We also observe that nothing changes if the quadratic part of the potential 



in (38) is omitted, so that also the purely quartic periodic chain with three 



particles is an integrable system. 



4 Explicit low-dimensional solutions 

From the expression ( |3^ ) of the numerical coefficients appearing in the general- 
ized force Fr in (^) it is easy to derive some selection rules already empirically 
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found in the literature. For instance, it is known (see e.g. ||^, |3^) that for a 
periodic FPU-/3 chain with an even number of oscillators the initial excitation 
of a set of modes all having even (resp. odd) indices only, cannot lead to the 
excitation of modes having odd (resp. even) indices. This property can be easily 
derived observing that initial conditions of these two kinds assign a special role, 
among all the coefficients C^jki i to those with one index of a given parity and all 
the others of the opposite parity. For example, if even modes only are initially 
excited, it follows from (36) that odd modes can be excited only through those 
Crjki with only one odd index; the converse applies to the initial excitation of 
odd modes only. Since the algebraic sums of indices of the A-symbols defining 
the coefficients in (^J) are, in both cases, always odd integers, they cannot be 
multiples of the even number N. Thus, all Crjki connecting (in the above sense) 
odd to even modes vanish. 

This example shows that (for even N), if we consider the set of all modes 
partitioned in the two subsets of the even and odd modes, this partition has 
the property that an initial excitation completely contained in one of the two 
subsets cannot propagate to the other. It is therefore natural to wonder about 
the existence of other partitions of the set of modes having the same property. 
Exploiting this idea we are able to find some explicit periodic and quasiperiodic 
solutions for the FPU-/3 system. Though their existence, for a finite number of 
oscillators N, is strongly dependent on the divisibility properties of N, once the 
solutions are obtained and expressed in terms of the original coordinates {qn} 
they remain valid for the infinite chain. 

Let us denote with Ai the set of all "internal" modes for an arbitrarily fixed 
number of particles N, i.e. the set of indices Ai = {1, . . . , A*" — 1}. The previous 
example leads us to introduce the following definition: a subset of modes A C Ai 
is defined to be "of type I" when 

Cjki =0 yi&M\A, yj,kJ&A. (40) 

It is clear that, if a subset A is of type I, every solution of ( ^6|) having Qi{0) = 
Qi{0) = 0, Vi S \ -A, is such that Qi{t) = OVt for the same indices i, which 
means that an initial excitation completely contained in A cannot propagate 
out of A. The remaining coordinates {Qj{t)}j^^ of the solution obey a reduced 
system of equations of motion analogous to (^6|) , with the only difference that 
the summation in the expression of the generalized force is restricted to indices 
j, k,l € A. So the normal coordinates and conjugated momenta with indices in 
A span an invariant subspace in the system's phase space, with dimension given 
by the double of the elements of A; the dynamics over it is generated by the 
reduced Hamiltonian 

HA{{Qj,PjheA) = IY.'^P-+^-Q1) + 

— ^ UJ^LOjLOk^^lCijklQiQjQkQl ■ (41) 
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The subsets of even and odd modes (for even A^) are just two examples 
of subsets of type I (each one with the additional non-generical peculiarity of 
having the complementary subset also of type I). It is easy to prove that these 
two subsets are type I independently of the chosen basis (p^). 

In the following subsections we discuss some low-dimensional type I subsets 
that we were able to identify, and their associated exact solutions. It must be 
remembered that we refer to the normal coordinates with 7 = 7r/4 in (^). 



4.1 One-mode solutions 

Let us begin with the case of the simplest subsets A: those consisting of only 
one mode, denoted by the index e. Then, the condition ( |40| ) for A = {e} to be 
of type I is 

Cneee =0 Vn G with 71 =^ e . (42) 

If this happens for some e, the corresponding reduced Hamiltonian (^) repre- 
sents a one degree of freedom (and thus integrable) system, described by the 
single coordinate Qe- 

Writing the required coefficients in terms of the A-symbols of ( ^2|) 



Cneee — — ^n+3e + 3A„_e , (43) 

we find that for each fixed e E Ai, except for 



N N N 2N 3N 

T' y' Y' ~3~' ~ 



(44) 



there is in M always one (and only one) n ^ e such that Cneee 7^ 0; thus 
condition ( ^ ) is not satisfied and {e} is not of type I. Given a mode e G M 
different from those listed in (^), the unique n €z A4 different from e such 
that Cneee ^ is in fact given by the unique integer in A4 congruent to — 3e 
modulo N. Such a mode can be called the mode n — n{e) "directly forced" by 
mode e, since in the generalized force Fn in ( ^ ) a forcing term proportional 
to CneeeQe is present. This does not mean that this mode will be the fastest 
growing mode, when mode e is the only one initially excited (since, of course, 
resonance relations must be examined), but its role is important because it acts 
as a trigger to excite other modes Qj (possibly more unstable but not directly 
forced by mode e) through terms of the kind CjneeQnQt in their equations of 
motion (|3^). 

On the contrary, for each of the modes listed in (^ ) (of course when some of 
them exists, i.e. when N has the divisibility property required for the considered 
e in (^4|) to be an integer) property (^) is verified; thus these are the only values 
of e such that {e} is a type I one-mode subset. If one of the modes in (^ ) is 
the only mode initially excited, it remains excited without transferring energy 
to any other mode, and the only modes that possess this property are those 
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in (Q). In this case, the equation of motion for the excited mode amphtude Qe 
is 

Qe - -^iQe - , (45) 

where Ceeee is always positive. The general solution of (^), having as free 
parameters the amplitude A and the time origin Iq, is 

Q,{t)=Acn[nS~to),k] , (46) 

where He and the modulus k of Jacobi elliptic cosine functiorj^ both depend on 
A: 

rie = IVeVrTS^ (47) 

with Se = /3wgCeeee/(2iV). We remind that this solution is periodic, with am- 
plitude dependent period Tg given in terms of the complete elliptic integral of 
the first kind K{k) by 

4:K(k) , 

T. . ^ . (49) 

Since A is one-to-one related to the control parameter involving the energy 
density, all parameters of the solution ( p6| ) can be rewritten as functions of e 
instead of A. For the period we have 

^ ^ (50) 



where k = k{e). The period decreases as e is increased, while in the limit e ^ 
it reduces to the harmonic period Tg — > iir/uje. 

All one-mode solutions correspond to stationary waves of the form 

g„(t) = Qe{t)w^:^ (51) 

where the vectors w^'^) = w^^^(7r/4) are given in formula (|l]). The spatial 
pattern is therefore the same of the corresponding linear mode, but these are 
nonlinear modes since the time dependence corresponds to a nonlinear oscilla- 
tion, which reduces to the usual harmonic one as A goes to zero. 

Let us first discuss the lattice pattern of the solution corresponding to e = 
N/2, i.e. its expression as a function of the old set of coordinates: 

qn{t) = ^yd^QN/2{t) , (52) 
V iV 



*For standard notation and properties of elliptic functions and integrals we refer to ]43| [i^. 
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where QN/2{t) is given by Eqs. (^), and (||), with e = cjjv/2 = 2, 

<^Ar/2 = 4/3/-/V. The lattice pattern has adjacent particles with opposite phases, 
being the same of the zone-boundary phonon mode. Actually, this solution is 
present for any potential of the form 



N 
n=l 

with even iV on a periodic lattice (or on an infinite lattice). In fact, after the 
substitution of the ansatz solution qn = (— l)"r(t), the equations of motion 
reduce to the single equation for r{t) 

r = -2 (2r) , (54) 

dr 

where ^evenif) = (^'('') + ^(~^))/2 is the even part of the $ function. For 
instance, a solution of this kind is known for the integrable case of the Toda 
lattice and corresponds to the standing cnoidal wave with two lattice spacings 
wavelength p5| |. 

The cases e — N/A and e = "iN/A are treated in a similar way. The lattice 
pattern is 



1 



qn{t) = ^Qe{t) ±sin(^— j +cos[— j (55) 



N 

where Qe{t) is given by Eqs. (^, and (^) with uj^/i = i^zn/a = V^i 
= ^■iN/A = 4/3 /A^, and the ' + ' and ' — ' signs correspond to e = and 
e = 2>N/A respectively. Let us observe that the solution e = 3A^/4 is obtained 
from that with e = N/A by shifting the lattice pattern by one lattice spacing; 
another shift produces the same solution apart from a temporal phase. This 
property is a consequence of the invariance of the model under translations of 
an integer number of lattice spacings. 

Because of an additional peculiarity, the single mode solutions e = N / 3 and 
e = 2 A/3 will be treated in the following as particular cases of the two- mode 
manifolds. 



4.2 Two-mode manifolds 

We do not try here to identify all multi-mode, or even two-mode, type I subsets 
for generic A^, although some of them can be found by simple inspection of a 
computed table of non-zero Ciju coefficientfj^ for small A^. Instead we confine 
ourselves to the search of two-mode type I subsets of the kind {i. A — i}, since in 
this case we get complete results leading to periodic and quasiperiodic explicit 
solutions. 

^For example, it can be seen that {1,5} and {3,7} for = 8, or {1, 5, 9}, {2, 6, 10} and 
{3, 7, 11} for Af = 12 are all of type I. 
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A subset {«, N — i} is of type I when (see condition (pO|)) 

Cnjki^O yneM\{i,N-i}, yj,k,l£{i,N~i} . (56) 

Thus, in particular, a necessary condition is 

Cniii = yn e M with n ^ i, N ~ i . (57) 

This can happen, a priori, only in two cases: a) {i} is of type I, or b) {i} is 
not type I but the mode n{i), directly excited by z, is coincident with N — i. 
But it is easy to see, from the explicit expression of n{i), that case b) can never 
happen: then if i is such that {i, N — i} is of type I, it must be one of the modes 
in (^ ) . The only subsets of the kind requested that can be type I are therefore 
£i = {N/4:,SN/A} and £2 = {N/3,2N/3}. To check that they are indeed, it is 
useful to remark that 



CN-i, N-j, N-k, N-l — Cijkl ■ (58) 

Then, since i verifies (42) with e — i (and also N ~ i does), from (58) we need 



only to check that 

Cn,^,^^N^^^O yn^M\{l,N -l} (59) 



for condition ( p6D to be satisfied. 

Let us first consider the set £i . Using formula ( |3^ ) we have 

Cn, N/A, N/4, 3N/4 — —\i+ZN/A + "^^n-^N/A + ^n+N/A j (60) 

and each A-symbol vanishes for n ^ £\, thus satisfying (|59|). This proves that 
the set £\ is of type I, i.e. the excitation of the the pair {iV/4,3iV/4} does not 
propagate to other modes. Moreover, Hamiltonian (^) with A = £\, describing 
the reduced system containing only this pair of modes, is separable because the 
coefficients (|60| ) are zero also for n <E £i, giving 



2 (^^N/A + -P3JV/4) + Q%/A + QIn/A + ^ (Qn/A + QsN/a) (61) 



Being separable, this Hamiltonian is integrable. Thus, the motion in the sub- 
space of the modes {N/i, 3iV/4} is simply given by the cartesian product of the 
one-mode solutions already introduced in the previous subsection. However, it 
is interesting to observe that this two-mode invariant manifold foliated in two- 
dimensional tori exists in the 2A''-dimensional phase space of the system, with 
periodic or quasiperiodic trajectories whose winding number depends on the 
non- linear periods of the two modes (see formula (^)). Therefore we are able to 
obtain quasiperiodic solutions with two frequencies exploiting the fact that an- 
harmonicity, causing the nonlinear dependence on amplitude of the period ([49|), 
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"removes" the harmonic frequency degeneracy of the pair 3A^/4}. 
general solution (see formula (E^)) is 



Quit) 



The 



(62) 



-j=Apf/4 cn [fJAr/4(Ajv/4) t + (j), kiAN/i)] 
V -/V 



sm I —n 



'TT 

COS I —n 



+ 



^=^3Ar/4 cn [Vt'iN/iiAzN/i) t + tp, ^(^3^/4)] 

V A" 



'TT 

sm I —n 



where we have emphasized the dependence of the various parameters on the two 
independent arbitrary amplitudes and A^]^/^. It is interesting to observe 

that 

qn+2 = -g™ , (63) 



which explains why a mode with a wavelength of four lattice sites corresponds 
to a solution with two degrees of freedom. When the amplitudes of the two 
modes in (^) are equal (we omit in this case the moduli in the arguments of 
the cn functions since they are equal) one gets periodic solutions. In particular, 
for the two choices of the phase difference ■0 — (/) = and — = 2K, stationary 
wave solutions are obtained. They correspond to a lattice pattern having one 
node every two lattice sites, thus giving also a solution valid for fixed end chains, 
after an appropriate choice of the lattice origin. It should also be observed that 
these latter would be one-mode solutions in the normal coordinates relative to 
the basis (|l^ ) with 7 = 0. 

If instead we choose ip — (j) = ±K [still with equal amplitudes, denoted by A, 
and with K = K{k{A))] the corresponding periodic solutions have the property 



qn+i{t ± At) = qn{t) , for At = K/nN/4{A) ,yn , 



(64) 



where the sign in the argument in the l.h.s. is the same of {tp — (f). We con- 
sider (|6^ ) as the defining property of a travelling wave on a lattice, propagating 
respectively to the right (' -I- ' sign) or to the left (' — ' sign) with velocity 1/ At 
(the lattice spacing being the unit of length) . We remark that the velocity de- 
pends on the amplitude, as should be expected for a nonlinear wave. Actually, 
observing that, for integer n, the following identities hold 

sin ^— cn(a;) = - {cn[iir(n — 1) + 2;] + cn[Jf(n — 1) — x]} , (65) 



cos yi^^j cn(a;) = 
and defining the function 



— {cT\{Kn + x) + cn(_R'n — a;)} 



(66) 



g[x, k) = cn(x — K{k), k) + cn{x, k) , 

it is straightforward to rewrite the general solution (|6^), for excitations in the 
£1 subspace, in a form which reveals its progressive and regressive travelling 
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wave content, in analogy with the harmonic case: 



qn{t) = (67) 
[Aag{Kan + riat + (j), ka) + Abg{Kb{n + 1) + ribt + iP, h)] + 

[Aag{Kan - nj - (f), ka) + Abg{Kb{n + 1) - J^bt - V, h)] ■ 



2VN 
1 



2VN 

In ( |67| ) we have denoted the quantities relative to the modes N/A and 3N/4: by 
the indices 'a' and '6' respectively. The particular cases above, with Aa = At, = 
A and ip — (p — ztK are 

qr,{t)^^g{KnTnN/4{A)tT4') , (68) 
V Jy 



where we have made explicit the travelling wave nature of these solutions. 
Turning now to the other set £2, we have 

Cn, N/3, N/3, 2N/3 = —^n+iN/S + 2A„_2Af/3 + , (69) 

and again we find that all A's vanish for n ^ £2, thus proving that also £2 
is of type I. This time there is a non-zero case among the different kinds of 
coefficients coupling the two modes of this subset, e.g. ^v/s. 2Ar/3, 2Ar/3 = Ij 

and the reduced Hamiltonian is 

^£2 = 2 i^^N/3 + ^2^7V/3) + 2 ('3w/3 + Q2N/3^ {Q%/3 + QIn/S^ ■ C''^) 

Let us first remark that this Hamiltonian reduces to the one of formula ( p8[ ) 
for iV = 3 apart from the center of mass motion term, which has already 
been eliminated. Being the potential in (|7^) invariant under rotation in the 
{Qn/3t Q2N/3) plane, again, as for the = 3 case, the "angular momentum" 
{QN/3P2N/3 — Q2N/3PN/3) is conserved. Thus, having two constants of mo- 
tion in involution, Hamiltonian (|70| ) is integrable. We have again an invariant 
manifold foliated by two-dimensional tori, but the explicit general solution for 
Hamiltonian (|70|), although possible in principle, is less trivial in this case. How- 
ever, we can easily find some lower dimensional solutions. In fact, the presence 
of a central potential implies that there exist solutions whose trajectory on the 
{Qn/3tQ2N/3) plane is an oscillation lying on an arbitrary chosen straight line 
crossing the origin. Hence, this implies that not only the modes N /i and 2N /Z 
are particular one-mode solutions, as we have already anticipated, but that 
there are analogous solutions corresponding to all possible choices (rotations) 
of the basis ( [l8| ) ; this property is not present for the set £1 . Each one of these 
solutions with zero "angular momentum" would be regarded as one-mode if we 
had chosen a linear basis with the convenient value of 7 in formula (^). In the 
original lattice coordinates they are stationary waves, with a pattern given 
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by an arbitrary continuous translation of the stationary sine wave of wavelength 
three. All these one-mode solutions satisfy the same differential equation (p5|), 
with uje — ujn/3 = \/3 and Ceeee = CAr/a, jv/3, w/s. w/s = 3. The lattice pattern 
is 

gnW = y^^cn [f7^/3(t-io), A:] sin + ' ^'^^^ 

where 7 is an arbitrary real parameter; fl^/s and k are given by ( ^ ) and ( ^Sf ) 
with de — Sn/3 — 9/3/(2Ar). It is interesting to observe that for a generic value 
of 7 this solution has no nodes (still oscillators) at lattice sites. The only case in 
which nodes appear is 7 = (the cases 7 = 7r/3 and 7 = 27r/3 can be reduced 
to this one after a translation of a lattice unit), and it is clearly present also for 
fixed ends oscillator chains. 

Also for the set £2 travelling wave solutions exist. They correspond to closed 
orbits for the central potential in Hamiltonian ([70|). A particularly simple case 
is the one of closed circular orbits: 

Q JV/3 (i ) = ^ COS [r2 jv/3 (i - io )] 



'^2N/3 

The lattice pattern is 



±Asin[r!Ar/3(t-to)] ■ ^ ' 



where 4> is an arbitrary phase. This solution is of the form of the linear 
travelling wave, but has amplitude dependent frequency and velocity, since 

Finally, we further stress that, though for the FPU-/3 system with finite N 
each of these solutions exists only for suitable (but still infinitely many) values 
of TV, their expressions in terms of the original {qn} coordinates remain valid 
as spatially periodic solutions for the infinite chain. The "out-of-phase" oscilla- 
tion ( p2[ ) and the one-mode iV/4 solution with 7 = (given by our solution ( |62[ ) 
in the particular case = A3JV/4, ip = (j)) appeared first in |4^, where 

they were obtained directly in particle displacement coordinates by a suitable 
"ansatz" . A statement equivalent, in our notation, to the fact that an excita- 
tion in the subsets £1 or £2 cannot propagate out of these sets, was given in 
unpublished work by ShoU [El) . 



5 Stability of periodic orbits 

In the previous Section we have shown the existence, in the phase space of the 
FPU-/? system, of various kinds of initial conditions giving rise to trajectories 
lying on manifolds of smaller dimension than that of the constant energy shell. 
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We have considered in detail some simple cases leading to reduced integrable sys- 
tems, thus obtaining some families of explicit solutions. In our opinion these are 
interesting in their own, since explicit analytical solutions, are not common at 
all for nonlinear non-integrable chains of oscillators. However we must remember 
that all these results hold for particular classes of initial conditions. Therefore 
an important problem is that of the stability of these solutions or, more gener- 
ally, of these invariant subspaces. If one of them turns out to be stable, then its 
qualitative features are present for a non-zero measure neighbourhood of initial 
conditions in phase space, thus acquiring physical significance. 

From another point of view, it is reasonable to expect some relation between 
loss of stability of periodic orbits (or more generally of invariant subspaces) and 
the onset of chaos in the system. This is suggested by the generic behaviour 
of near-integrable Hamiltonian systems where, as the perturbation parameter 
is increased, the "chaotic sea" expands in phase space, forcing regular regions 
to shrink. Budinsky and Bountis conjectured that the destabilization of 
known periodic orbits of an Hamiltonian system, although being only a local 
property, could give an upper bound estimate for the onset of large-scale chaos. 
This suggestion was perhaps not fully justified since, at that time, only few one- 
mode solutions (e.g. e = N/2) were known, and moreover the behaviour of very 
low-dimensional systems (e.g. Henon-Heiles) had been studied. It assumes more 
relevance now, since we know that several invariant manifolds exists. However, 
the fact that periodic and quasiperiodic solutions lose stability or invariant 
submanifolds become unstable does not necessarily imply large-scale chaos. We 
will show that this loss of stability leads to the formation of stochastic layers, but 
their thickness and extension in phase space remain small at small e. In modal 
space, the instability of all these newly discovered submanifolds does not lead to 
energy equipartition if e is sufficiently small, but to exponentially localized mode 
energy spectra with irregular, chaotic, time dependencies of mode energies. 



5.1 Bifurcation and transition to chaos in a reduced sys- 
tem 

It is easy to verify that the set T = {N/A, N/2, 3N/A} is of type I. Again, 
one has to verify condition ( pO| ) for A = T. Apart from permutations of the 
indices, there are ten distinct kinds of coefficients to consider, corresponding 
to the number of triplets that can be formed with the elements of J- without 
regard to order and with allowed r epet ition of any element. Five of them belong 



to cases already treated in Sees. and |4.2| . Examining also the remaining 
cases, we find that the only non-zero coefficients having three indices in 
and the remaining index unrestricted (in M) are, apart from permutations, 

CN/i,N/i,N/2,N/2 — ^^3^/4,3^/4,^/2,^/2 — 2 and, of COUrSC, the three Ceeee 

with e ^ T . This allows to conclude that T is of type I, and at the same time 
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leads to write the corresponding reduced Hamiltonian: 

Hj^ = - + Pn/2 + -^3^/4) + Qn/4 + '^Q'n/2 + QIn/4 



N 



Either from the corresponding equations of motion, or from what has been found 
above concerning the coefficients, we observe that there exist also two new two- 
mode invariant subspaces: that of the pair {N/A, N/2} and that of {3N/4, N/2}. 



Qn/4: — ^'^Qn/4 — 8-^Q'^/4 — 24-^(5^/4(5^^2 

Qn/2 = ^4:Qn/2 ~ 16-^(5^/2 ~ 24-^(5 7V/2QAr/4 



At variance with cases discussed in Sees. 4.l|^4.2|, all these new invariant sets 
give rise to seemingly non-integrable Hamiltonians. The case of a two-mode 
invariant subset is particularly interesting because, for such a reduced system, 
the qualitative features of the dynamics can be visualized on a two-dimensional 
Poincare section. 

The equations of motion corresponding to the set {N/A, N/2} are: 

(75) 

Those relative to the other set {3N/4, N/2} are identical, substituting Qn/4 
Qs,N/4- Rescaling the coordinates by (3/N one obtains equations independent 
of the ratio (3/N. Thus an orbit of the system ( |75| ) is mapped to one of the same 
system with /3 /N = 1 . This means that the system (^5|) remains meaningful in 
the termodynamic limit. 

We want, in particular, to investigate the stability of the "out-of-phase" 
nonlinear oscillation (|5^), i.e. the one- mode periodic solution ( ^6|) having e = 
N/2. Therefore we choose the {Qn/4, Pn/4) plane as the surface of section; in 
fact on this plane that solution is represented by a fixed point for the Poincare 
first return map, located at {Qn/4i Pn/4) = (0,0). The one-mode solution with 
e = N/4: lies, instead, on the surface of section and can be drawn simply as the 
curve bounding the energetically accessible region. In Figs. 1, 2 and 3 we show 
numerically generated Poincare sections for the system ([ts]), corresponding to 
three different increasing values of the control parameter e in (|3|). For each 
e value, the intersections originated by several different initial conditions are 
shown, together with the bounding curve corresponding to mode e = iV/4. 

We observe that the fixed point at the origin is stable at low values of e. 
When the energy density is increased, a bifurcation occurs, giving rise to a 
pair of stable fixed points, while the original one loses stability. A small closed 
invariant curve which surrounds only one of them is plotted in Fig. 2 to make 
clear the presence of two distinct fixed points (i.e. two different bifurcated orbits 
for the flow), so that the bifurcation is not period-doubling. To confirm this we 
have also numerically computed the eigenvalues (Ai,A2) of the Poincare map 
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linearized around the origin, obtaining a more precise estimation of the critical 
e value. As e is increased the eigenvalues, lying on the unit circle in the complex 
plane, move on it until they collide at Ai = A2 = +1 for e = 1.127 and then 
move apart on the real axis. After the fixed point at the origin turns unstable, a 
stochastic layer is originated by homoclinic intersections, too thin to be seen in 
Fig. 2 but already very thick in Fig. 3. Therefore when the stationary wave ( ^2[ ) 
has high enough amplitude and is perturbed only by a small component along 
mode N/4:, the interaction between the two modes leads to chaotic behaviour. 
However even at the e value of Fig. 3 we do not observe transport over the full 
phase space which is signature of large-scale chaos (which, of course, will set in at 
larger values of e) . This shows that the loss of stability of periodic orbits induces 
the appearence of stochastic layers but the large-scale chaos phenomenon is not 
a direct consequence of it, contrary to the expectations of Ref. 

In the next Subsection we study the stability problem for the e = N/2 mode 
in the full phase space. 



5.2 Large N stability 



The problem of stability of a known solution Q{t) under generic small per- 
turbations in the initial data, is tackled by studying the variational equations, 
obtained by linearizing system (36) around Q{t). Let x{t) = Q{t) ~ Q{t) be 
the small separation vector between a solution Q{t), initially close to Q, and 
the reference solution Q{t) itself. Then x obey, in linear approximation, the 
following system of equations: 



N-l 



Brs{t)Xs, 



r = 1, 



,7V- 1 



(76) 



where 



Brs{t) = 



dFr 



dQs 



Q=Q(t) 



(77) 



is the Jacobian matrix of F evaluated on the reference solution. We limit 
ourselves to study here the stability of the periodic single-mode {e} solutions 
Qj{t) = SjeQe{t), with Qe{t) givcu by (^. In this case system (^ reduce to 



3/3 

r = -UJ^Xr - —UJrUjlQl{t) ^ W,C<;^ 



Af-1 



2A^ 



r = 1, 



,N-l 



(78) 



where the coupling matrix, which depends on the chosen mode e in (Q), has 
elements Cr^s = Crsee- It turns out that c'fs is diagonal only for e = A^/2, 
with ci^^'^^ — 2drs- The study of the stability of the e — N/2 nonlinear mode 
is therefore simpler, because the different components of the perturbation in 
modal space are all decoupled and can be studied separately, each one being 
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a problem of parametric excitation described by an equation of the Hill type. 
This permits approximate semi-analytic estimates also for what concerns the 
dependence of stability properties on the number TV of oscillators. Accordingly, 
we restrict here to the case e = N/2. The study of other one-mode cases (e.g. 
e = would imply the study of the stability for a system of coupled linear 

equations with periodic coefficients, a more complex problem of many degrees 
of freedom parametric excitation. This could be done by producing numerically 
the appropriate transfer matrix, but such a procedure should be repeated for 
each fixed value of N and the large N result should be extrapolated. 



For e = iV/2, after using (46), Eqs. ( [ZSj ) reduce to a set of decoupled Lame 
equations Q 



1 

TV - 1 , (79) 



dA 

l + 12^cn2(17jv/2i,fc) 



with f2Ar/2 and k depending on A as given in (1^) and (ll). The physical con- 
trol parameter, proportional to the energy density, is e, which is in a one-to-one 
relation with A. However the analysis is simpler if we choose the modulus fc as 
a control parameter. In fact, as e is varied from zero to -foo, the modulus fc in- 
creases monotonically in the interval [0,1/ \/2 ) , the limit e — > +oo corresponding 
to fc ^ \/V2. The relation e{k) is 

'l + T^Ij) ■ (80) 



1 - 2fc2 y 1 - 2A:2 

Expressing and PA^ /N as functions of fc, and after rescaling time at each 

fixed k (which does not alter stability properties), each of the equations in (Foh 



can be put in the standard Jacobian form of Lame equation |47, ^ 

y" -t- [a - v{v -t- l)/e2sn2 (m, fc)] y = , (81) 

where the prime superscript denotes differentiation with respect to the new 
time u = ^N/2t and y stands for the generic variable Xr- For the perturbation 
component Xr in the direction of mode number r, the parameters in (pll) are 



a^{l + Ak^)uj^j4; + 1) = 3lu^. /2 . (82) 

We observe that dependence on the mode number r and on the number N of 
oscillators comes only through the ratio r/N appearing in ujr (see eq. (|l3|)). 
This allows to discuss in an unified way also the dependence on the system size 
N. In fact, we can fix e and then consider two systems of different sizes A^i and 
the linearized equations for a mode ri in system A^i and for a mode r2 in 
system N2 are the same if ri/Ni — r2/A^2- 

Resuming, we discuss the stability of the zero solution of ( ^ ) as a function 
of two parameters (p, k) with 0<p<l, 0<fc< l/\/2, when {a, v) are given 

by 

a = p(l + 4fc2); y{v + l) = Qp. (83) 
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A given pair (p, k) describes at the same time different mode numbers r in 
systems with different N, all having the same r/N value given by the relation 
p = sm^{TTr/N), with the same e for the unperturbed solution e = N/2. We first 
consider p as a continuous parameter, which corresponds to the limit N — > oo, 
and then we discuss what happens when p can take only a finite set of values, 
as it is for any finite N . 

For integer v many rigorous results are known for eq. (|8^) |4^, and it 
is in particular known [^9| that there are only v + 1 instability regions in the 
(a, k) plane. However, in our case only the integer values 1^=1 and i' = 2 
are possible, corresponding to p = 1/3 and p = 1, respectively. The stability 
charts in the (a, k) plane for these two integer v cases of eq. (|8^) have been 
explicitly constructed in [Q. The second case p = 1, is of limited interest to 
us since it corresponds to r = N/2, i.e. to perturbations along the invariant 
subspace e = N/2; the results in just confirm that there exists a periodic 
solution for xpj^2 ^lII k values, as it must a priori be, since it corresponds to 
a perturbation along the periodic orbit QN/2it) under study. More interesting 
is the case p = 1/3; the representative point in the {a,k) parameter plane 
of dsT] ) with V — 1 moves on the curve a = (1 + 4fc^)/3 as k is varied, and 
we find (see Fig. 4 of Ref. [Q) that it lies always in the first stable region 
for < fc < l/\/2. It touches the boundary curve of the instability tongue 
only for k = 1/2. Therefore the mode corresponding to p = 1/3 is stable (as a 
perturbation to e = N /2) at any finite energy density e, and reaches the stability 
border only asymptotically for e — > oo. 

To solve the stability problem in the generic v case we approximately reduce 
eq. ( |8l| ) to Mathieu equation, whose stability diagram is well known Q, 
In fact, in the parameter region of our interest < k^ < 1/2, the function sn^ 
appearing in (|8^) is well approximated by its first order Fourier development. 
This latter is given by 

fc^sn^ M, fc^l — > ^cos , 84) 

^ ' ' K \K) ^^siBh{m:K' /K) \ K ) ' ^ ' 

where K and E are the complete elliptic integrals of first and second kind, 
respectively, with modulus fc, while K'{k) = K{\/1 — k"^ ). The Fourier coeffi- 
cients are exponentially decreasing with order n, the more rapidly the smaller 
k is. In the worst case to consider, i.e. k 1/V2 the ratio between the n = 2 
and n — 1 coefficients is « 0.086. Furthermore, we are interested in the modes 
which first go unstable, i.e. which possibly destabilize at small fc, and for them 
the error is much smaller. Therefore, we truncate ( |8^ ) keeping only the n = 1 
term and obtain, after a further time rescaling r = ■ku/2K, the canonical form 
of Mathieu equation: 

^ + [a-2gcos(2T)]2/ = , (85) 
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where 



(- + «^ + 4)^ '-SOT)' 

Both parameters a and q depend on fc^ and p, therefore for each fixed p (which 
corresponds to fix a ratio r/N), changing (and correspondingly the energy 
density e) amounts to trace a curve in the (g, a) parameter plane. A stability 
transition happens if this curve intersects characteristic curves separating stable 
from unstable regions for the Mathieu equation. In our case intersections can 
occur only with the characteristic curve 61(g) (we use the standard notation, 
see e.g. @) bound ing from below the first tongue of instability. In Fig. 4 we 
show the (g, a) curves (^) at fixed p and varying k for positive^ g, and those 
at fixed k and varying p, together with the characteristic curves 61(g) and ai(g) 
bounding the first instability tongue. This mesh provides a reference frame to 
locate the values of p and k corresponding to stable or unstable cases. We do 
not trace the p = 1 curve which superposes with ai(g), showing only a slight 
deviation for > 0.4. This superposition is expected because the truncation 
to the first Fourier coefficients in (|^ only slightly affect the exact solution of 
the Lame equation in the v — \ case. We see that curves with fixed p enter 
the instability region, for some critical value kc{p) of /c, only for p larger than 
a value such that the corresponding curve touches 61 (g) only asymptotically for 
k 00. We numerically find such a value to be p ~ 0.33, which therefore 
must be identified with the value p = 1/3 having rigorously this property for 
the exact Lame equation, as seen before. This confirms the accuracy of the 
Mathieu approximation. Also, we find that the p — 1/2 {r = N/4) curve 
intersects 61 for e ~ 1.13, in excellent agreement with the value e — 1.127 for 
the instability threshold in the r — N/A direction, obtained from the Poincare 



map in Subsection 5.1 



We have computed, for several values of p > 1/3, the critical kdp) value 
at which the corresponding constant p curve intersects the 61 line. Returning 
to the physical control parameter e by ( pO[ ) we finally construct the marginal 
stability curve edp) in Fig. 5. Now we remember that, for a given number N 
of oscillators, only the discrete set of p values 

p = sin^{Trr/N) , r = l,...,7V-l (87) 

must be considered. From Fig. 5 we see that for each mode having p > 1/3 there 
is a threshold value edp) ^ for instability. Above edp) the e = N/2 nonlinear 
mode develops an instability causing growth of the mode corresponding to p 
through parametric resonance. On the contrary, modes with p < 1/3 (i.e. r/N < 
0.196) are always stable in the linear approximation for any energy density of 

^We used positive q values for the curves ( |86[ ) since the Mathieu stability chart is symmetric 
about the a-axis. 
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mode e = N /2, so that (nongeneric) perturbations of the latter involving only 
these modes never lead to instability. These modes, as well as modes with 
p > 1/3 when e < ec(/c), can grow only if they are triggered by the interaction 
with other modes that are unstable. This kind of interaction is neglected in the 
linear approximation around the e = N/2 mode, and comes into play only at 
later times, when unstable modes have grown and consequently the linearized 
theory is no longer valid. 

We can now describe what happens for a given system of N particles. First 
of all, for A'^ > 4 there are always modes with p > 1/3 so that the e = N/2 
mode can never be stable for all energy densities. As Cc is a decreasing function 
of p, the first modes to go unstable when e is increased from zero are always 
r = {N/2) - 1 and r = {N/2) + 1, which have p = cos^ {-k/N). Therefore for each 
(even) number N of particles there is a non-zero value e{N) = ec{cos'^ {tt / N)) 
below which the nonlinear mode e = N/2 is stable. In fact for e < e{N) all the 
discrete points (p, e) having abscissas (pj) are below the critical curve ec(p) in 
Fig. 5. Since Ec ^ for p ^ 1, the critical threshold e{N) for the stability of the 
nonlinear out-of-phase mode approaches zero as the number N of particles is 
increased without limit. Using power series expansions of the curve in (^) and of 
hi{q) around the point (g,a) = (0, 1), we find that i{N) = tt^ /{■iN^)+0{N-'^). 
This is in agreement with the result obtained in Ref. using infinite Hill 
determinants. However, we do not agree with their results in Table I, which 
seem to indicate that low modes become unstable, while we claim that all modes 
with p < 1/3 are always stable. 

Furthermore, as the p values (^) to be considered are increasingly dense 
as N is increased, for each e, no matter how small, there is a value N for the 
system size above which the e = N /2 mode is unstable. In this sense, we can say 
that the e = N /2 mode is unstable at any nonzero value of the energy density 
in the thermodynamic limit. However for very small e it mantains a metastable 
character, with a non-negligible lifetime, because the largest growth rate of 
unstable modes vanishes as e ^ 0. This can be understood from Fig. 4 since 
for fc ^ the unstable part of the constant k curve tends to the edge between 
the 6i and ai lines and thus corresponds to a smaller and smaller growth ratej^ 

To confirm our analysis we have numerically integrated the full nonlinear 
FPU model directly in normal coordinate space (Eqs. (^6|)). We have overcome 
the difficulty of the huge number of interaction terms by first calculating a table 
of nonzero interaction coefficient, whose number is only ~ N'^ , which is then 
stored in a single vector variable with a suitable label ordering. We excite 
mode e = N/2 with an energy -Ejv/2 corresponding to a chosen e and give to 
all other modes an equally distributed small amount of energy, e.g. such that 
Sj5^Ar/2 ~ 10^^^. Then we follow the time evolution of all the modes, 

observing in particular their harmonic energies. The latter are the physically 

^We point out that our previous scalings of the time variable do not change the order of 
magnitude of time scales at small e since t ^ 2t for e — > 0. 
*We thank J. Laskar for suggesting this trick. 
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relevant observables, and are observed to grow exponentially for unstable modes, 
while the corresponding normal coordinates oscillate with exponentially growing 
envelope. 

First of all we have verified the scaling suggested by (jS^) . To compare results 
with different N we mantain fixed the ratio J2j^N/2 so that changing 

the A'' value at fixed e let modes with the same r/N ratio start with the same 
normal coordinate amplitude. In Fig. 6 we show the modal energy spectrum at 
fixed time for various values of N, with the energy density associated to mode 
e = N/2 fixed at e = 0.1. The superposition of different spectra proves that only 
direct interaction between mode e = N/2 and other modes, through parametric 
excitation, is important in a first time evolution, when the only growing modes 
are those predicted to be unstable on the basis of the Mathieu equation stability 
diagram. Modes with r/N < 0.35 are predicted to be stable for e = 0.1 from the 
curve €c{p) shown in Fig. 5, in agreement with Fig. 6. We have observed that 
at later times stable modes begin to grow due to secondary interaction with the 
unstable modes grown meanwhile. 

To understand the e — > limit for the infinite chain we have numerically 
computed the behaviour of the fastest growth rate of the instability when e is 
decreased while N is increased. In fact, if e is decreased at fixed N, when it 
passes below e{N) no instability is detected. Therefore to measure instability 
growth rates as e decreases, we must progressively increase N. We have detected 
the fastest growing mode, fc, at each value of e, and measured its exponential 
growth rate A {E^, ~ 10^*). The dependence of A on e, for small e, can be 
numerically fitted with a phcnomenological law A — Cie'^ + C2 with d = 1.0 
Ci = 1.1 and \c2\ < 6 x lO^'' which linearly extrapolates to a value consistent 
with zero in the limit e — )■ 0. We have therefore an interpretation of what might 
happen in the thermodynamic limit: in this limit the mode e = N/2 is always 
unstable under generic small perturbations but the rate of instability goes to 
zero (a sort of marginal case). It is interesting to observe that interchanging 
of the two limits e — > and the thermodynamic limit N —^ oo gives different 
answers for the stability of the e = N/2 mode. Stability is obtained if e ^ is 
performed first, while performing first the N ^ oo limit we obtain instability. 

An interesting and open question is the fate of the instability, since our 
linearized theory describes the short time scale only. For large N, the energy 
density stability threshold e{N) is very small, and the e = N/2 mode is unstable 
for all relevant e values. However one should still distinguish different regimes. 
For e values large with respect to those of the instability, the system evolves 
towards equipartition, as it is in general known from studies on the threshold 
of strong stochasticity for the FPU-/? system, although very little attention has 
been paid to high mode initial excitations, apart from Ref. ||]. However, for 
small e, numerical results show that the (slow) growth of unstable modes sat- 
urates at later times leading to a characteristic asymptotic form of the energy 
spectrum where modal energies perform characteristic oscillations and recur- 
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rences (similar to the well known FPU recurrences present for long wavelength 
initial excitations). Furthermore the spectrum decays, roughly exponentially, for 
modes far from mode number N /2, and it might be possible the emergence of 
some coherent structure. In this connection we point out that the more unstable 
modes are those around the initially excited mode number, a characteristic fea- 
ture of modulational instability, well studied in continuum systems such as, e.g., 
water waves |51 . In the context of discrete systems different from FPU models, 
namely chains made of harmonically coupled particles subjected to a nonlinear 
on-site potential, modulational instability has been proposed |]5^ as the first 
step towards the creation of states with spatially localized energy. The possible 
existence of spatially localized oscillations, called intrinsic localized modes, in 
FPU-type chains has been the object of much work, starting from gener- 
ally unrelated with studies on equipartition and energy thresholds in the same 
systems. Let us finally mention that in a recent interesting paper by Sandusky 
and Page a connection is investigated between the existence of intrinsic 
localized modes and the instability of the "out-of-phase" mode in several kinds 
of one-dimensional homogeneous chains with anharmonic intersite coupling. 



6 Conclusions 

There has been recently a revival of interest in conservative coupled oscillator 
systems (see e.g. |5^, ^). In this paper we have concentrated our attention 
on the prototype of these systems, the Fermi-Pasta-Ulam oscillator chain Q. 
Although this system has been studied for a long time, still it continues to 
be the source of many interesting physical problems. The first Section of this 
paper is devoted to a comprehensive review of old and new results (a recent 
historical reconstruction can be found in Then, we turn to linear mode 

space analysis. This was first performed in the basic paper by Izrailev and 
Chirikov [||, although the scope was there different. We have here derived 
a rather compact expression for mode coupling coefficients (see formulas ( ^2|) 
and (34)), which leads to an efhcient writing of the Hamiltonian in mode space 
(see (33) and (|35|)). After an "intermezzo" devoted to the iV = 3 system, where 
we rederive a well known integrability result by a one line proof, we pass to 
the central issue of the paper: the existence of low-dimensional exact soutions 
and invariant submanifolds. The method for proving their existence relies on a 
detailed analysis of mode interaction coefficients and of selection rules which pre- 
vent direct interactions among particular subsets of modes. We have thus found 
one-mode and two-mode integrable Hamiltonians (the latter because separable 
or rotation invariant) corresponding to new standing and travelling nonlinear 
wave solutions [see (^ , (|5|) , (|2|) , (|6^) , (|8|),(|7l|),(^] of the FPU system, which 
hold true in the thermodynamic limit. Three-mode invariant submanifolds are 
moreover present (see (fz^)), for which a Poincare section study of two-mode as- 
sociated subspaces reveals a bifurcation as the control parameter e, proportional 
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to energy per oscillator, is varied. This bifurcation is induced by parametric 
perturbation and leads to the formation of stochastic layers, corresponding to 
intermode energy transfer. This effect suggests a physical interpretation of these 
bifurcations, which could be responsible for opening channels of energy transfer 
among modes. It should however be reminded that all submanifolds so far found 
correspond to high k modes. The mechanisms controlling energy transfer for 



low k modes seem to be completely different (see Refs. |33| [24|, pO], p6[) 



Moreover, the opening of a transfer channel does not necessarily mean that 
energy flow is completed until energy equipartition is reached; this can happen 
only if large-scale chaos is present, leading to the diffusion of the orbits over the 
whole phase space. 

A complete stability analysis is possible for mode e = N/2, after reducing 
variational equations ( [zs] ) to a single Lame equation (|8l|) and a proper interpre- 
tation of control parameters related to mode number, system size N and energy 
density. For finite TV we prove the existence of a low energy range of stability. 
However this range shrinks to zero as TV ^ oo but, correspondingly, the rate of 
instability goes to zero as e 0, thus proving the presence of a sort of marginal 
stability of mode e = N/2 in the thermodynamic limit at low energy density. 
This is a hint to the belief that this solution might be physically relevant, and 
moreover encourages to look for analogous stability proofs for other one-mode 
and even two-mode solutions, whose existence was proven in this paper. 

A more difhcult matter would be the proof of stability of invariant subman- 
ifolds. The odd and even invariant submanifolds mentioned at the beginning 
of Section ^ also display an instability threshold. This phenomenon has been 
studied in connection with an instability of cnoidal waves in the modified KdV 
equation (mKdV) by DriscoU and O'Neil the control parameter for this 
instability might be different from energy density, as claimed in Ref. . How- 
ever the relation with this problem is not clear because the mKdV describes only 
the long wavelength part of the mode spectrum. The study of the stability of 
multi-mode invariant submanifolds will be the subject of future investigations, 
together with a more complete computer assisted search and classification of 
them. 
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Figure captions 



Fig. 1. Poincare section for the system ( |75| ) on the {Qn/4:t Pn/4) plane, with 
f3/N — 0.4. The energy value E corresponds to e = (3E/N = 1.0. Each 
smooth curve corresponds to a different initial condition on the constant 
energy shell. The fixed point at the origin, representing the periodic os- 
cillation ( ^2|) on the lattice, is stable. 

Fig. 2. Poincare section for the system (|75| ) on the {Qn/4i Pn/a) plane, with 
P/N = 0.4. The energy value E corresponds to e = PE/N ^ 1.25. The 
fixed point at the origin has become unstable and a pair of new stable 
fixed points has bifurcated from it (see text for details). 

Fig. 3. Poincare section for the system ( |75| ) on the {Qn/a^Pn/a) plane, with 
j3/N = 0.4. The energy value E corresponds to e = (3E/N = 10.0. The 
fixed point at the origin, which is by now unstable, lies within a stochastic 
layer. 

Fig. 4. Stability chart of the Mathieu equation (|85|). The instability region lies 
between the continuous lines. bi{q) is the lower bound and ai{q) the upper 
one. The lines match at (0, 1). We draw in the figure also the (q, a) lines 
in formula (|86|), parametrized by p and k. The dashed lines correspond to 
a fixed p and a continuosly varying k from to 1/ ^/2\ they are traced at p 
intervals of 0.1 from p = Q io p = 0.9. The dotted lines correspond instead 
to a fixed k and a continuously varying p from zero to one; the curves are 
traced at k^ intervals of 0.05 from 0.05 to 0.5. Crossings of the lines in 
formula ( |86| ) with the hi{q) line correspond to transitions to instability. 

Fig. 5. Marginal stability curve ec(p) (full line). The vertical asymptote at 
jO = 1/3 (dotted) is also traced to mark the fact that all modes with 
p < 1/3 are stable for any e value. 

Fig. 6. Modal linear energies E^ vs. the ratio of the mode number r to the 
system size N at fixed time t — iQ. Different symbols correspond to 
different system sizes: (+) TV = 8; (x) iV = 16; (□) N = 32; (O) N = 64. 
The continuous curve is an interpolation to guide the eyes. The initial 
values of modal energies are 3.2 x 10^^^. 
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